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1.  Introduction 


A  basic  type  of  magnetic  guiding  structure  used  on  atom  chips  is  fonned  by  placing  a  long, 
straight  current-carrying  wire  in  an  externally  produced  uniform  magnetic  field  that  points  in  a 
direction  perpendicular  to  the  track  of  the  wire  (Fortagh,  2007).  In  this  configuration,  the  roughly 
circular  magnetic  field  of  the  wire  is  cancelled  along  a  line  that  parallels  the  wire,  creating  a  long 
straight  line  of  zero  magnetic  field.  Cold  magnetic  atoms  can  be  trapped  in  quantum  states  in  this 
type  of  magnetic  field  minimum  near  to  the  wire  and  be  guided  along  parallel  to  the  wire 
direction.  This  forms  a  type  of  magnetic  waveguide  for  atoms  that  can,  in  principle,  be  used  for 
constructing  atom  interferometers  and  other  devices  based  on  atom  guiding. 

In  an  atomic  waveguide  based  on  magnetic  interactions,  the  forces  on  different  spin  components 
have  opposite  signs.  Thus,  if  one  component  of  a  spinor  is  bound  near  the  guiding  center,  another 
component  of  the  spinor  is  repelled  from  the  center.  This  is  simply  an  example  of  the  Stem- 
Gerlach  effect  commonly  used  to  separate  particles  in  different  spin  states.  One  difference 
between  the  simple  Stern-Gerlach  experiment,  in  which  only  the  spin  degrees  of  freedom  are 
quantized  and  the  atom  waveguide,  is  that  the  spatial  degrees  of  freedom  must  be  quantized  as 
well  as  the  spin  in  the  atom  guide  problem. 

In  order  to  accurately  treat  the  complete  quantum  problem  both  inside  and  outside  the  guiding 
region,  a  power  series  solution  expanded  about  the  guide  center,  yet  valid  for  the  whole  real  axis, 
has  been  developed.  Using  this  series,  both  the  bound  and  unbound  components  of  the  full  spinor 
solution  can  be  treated  within  the  same  framework.  The  power  series  converges  quickly  at  very 
small  radii  but  converges  very  slowly  at  relatively  small  distances  from  the  guide  center. 
However,  this  limitation  is  seen  to  be  due  to  the  fixed  hardware  floating  point  resolution 
available  in  modern  computers.  The  power  series  solutions  can  be  readily  summed  over  a  wide 
range  of  the  radial  coordinate  by  making  use  of  software  floating  point  techniques  implemented 
in  arbitrary  precision  arithmetic  libraries.  This  technique  produces  high  precision  eigenvalues  as 
well  as  detailed  infonnation  about  the  behavior  of  the  spinor  wavefunction  at  both  large  and 
small  distances  from  the  origin,  making  it  possible  to  perfonn  detailed  studies  of  atom 
propagation  without  resorting  to  the  adiabatic  approximation  (Sukumar,  1997). 

Numerical  solutions  for  the  quantum  modes  and  eigenvalues  of  spin  1/2  and  spin  1  particles  in  a 
quadrupole  guide  have  been  published  elsewhere  (Hinds,  2000;  Potvliege,  2001;  Lesanovsky, 
2004;  Bill,  2006);  however,  in  all  of  these  works,  the  boundary  conditions  used  at  the  guide 
center  have  been  chosen  incorrectly  and  do  not  properly  treat  the  regular  singularity  that  exists 
there.  This  error  in  the  boundary  conditions  clearly  results  from  treating  the  two  second-order 
radial  equations  (our  equations  8)  as  though  they  uncouple  at  the  origin  before  trying  to  perform 
a  power  series  analysis  on  the  full  system  of  radial  equations.  The  power  series  analysis  can  be 
performed  without  first  decoupling  the  system  but  one  must  be  careful  not  to  assume  that 
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relations  like  pR_  (p) «  R+  ( p)  hold  as  p  — »  0 .  This  type  of  reasoning  seems  to  have  caused 

problems  with  the  boundary  conditions  used  in  these  earlier  works.  In  Blumel  (1991),  the 
Frobenius  method  was  used  on  uncoupled  equations  to  properly  handle  similar  types  of 
singularities  in  the  problem  of  a  neutron  in  the  magnetic  field  of  a  rectilinear  current,  a  different 
but  related  problem. 

Failure  to  handle  the  boundary  conditions  properly  leads  to  an  incorrect  description  of  important 
details  of  waveguide  properties  such  as  the  energy  spectrum  and  the  ground  state  of  the  system. 
In  addition,  a  whole  class  of  modes  goes  undetected  in  these  calculations  simply  because  the 
boundary  conditions  are  not  treated  carefully.  The  quantum  properties  of  an  atom  guide  cannot 
be  treated  consistently  using  an  inconsistent  basis  set.  Since  the  quantum  system  is  linear,  using 
incorrect  boundary  conditions  in  the  numerical  solution  leads  to  apparent  modes  that  are  actually 
superpositions  of  the  true  modes  of  the  system.  These  mixed  modes  cannot  generally  be 
stationary  states  of  the  system;  they  are  mixtures  of  nondegenerate  states.  Various  lifetime 
estimates  have  been  performed  by  several  authors  (Hinds,  2000;  Sukamar,  1997;  Potvliege, 

2001;  Lesanovsky,  2004;  Bill,  2006)  using  these  nonstationary  modes.  All  such  estimates  are 
questionable,  since  the  mixed  modes  used  for  those  calculations  would  be  time  dependent  simply 
due  to  the  normal  time  evolution  of  a  conservative  Hamiltonian  system  and  this  point  is  not  made 
clear  in  those  works.  Thus,  any  time-dependent  properties  calculated  using  these  nonstationary 
modes  might  actually  be  due  to  the  natural  unitary  evolution  of  the  conservative  system  and  not 
representative  of  an  inherent  limitation  of  magnetic  guiding.  This  needs  to  be  clarified  before 
attributing  a  lifetime  to  the  states. 

In  a  previous  work  (Golding,  2009),  series  solutions  to  the  radial  differential  equations  were  used 
simply  to  provide  accurate  initial  conditions  appropriate  for  the  regular  singularity  at  the  origin. 
However,  in  that  work,  the  initial  conditions  were  simply  used  to  initialize  various  hardware 
precision  numerical  integration  routines  that  are  routinely  used  for  solving  ordinary  differential 
equations  (ODE).  The  shooting  technique,  along  with  a  root  solving  approach,  was  used  to  solve 
for  the  eigenvalues  of  the  system  as  well  as  plot  the  radial  eigenfunctions.  The  differential 
equations  solved  are  stiff  and  it  is  difficult  to  judge  the  accuracy  of  the  eigenvalue  results 
obtained  using  the  shooting  technique.  This  is  a  well  known  problem  with  the  shooting  technique 
that  occurs  when  an  exponentially  growing  solution  is  present  along  with  an  exponentially 
decaying  solution.  When  the  value  of  the  decaying  solution  is  too  small,  round  off  errors  can 
cause  the  growing  exponential  solution  to  become  dominant,  limiting  the  resolution  of  the 
calculated  eigenvalues.  In  second  order  systems,  this  problem  is  often  handled  by  integrating  the 
system  in  both  directions  and  matching  the  two  solutions  in  a  region  where  both  are  judged  to  be 
accurate.  However,  in  the  fourth-order  system  studied  here,  there  are  two  oscillatory  solutions  in 
addition  to  the  two  exponential  solutions  that  contribute  significantly  at  large  distances  from  the 
center.  The  oscillatory  solutions  cause  practical  difficulties  in  establishing  satisfactory  boundary 
conditions  for  an  inward  integration  of  the  system  since  one  cannot  simply  assume  that  the 
solution  simply  decays  to  zero  far  from  the  guide  center. 
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In  this  work,  the  solutions  to  the  quadrupole  guide  problem  are  formulated  using  a  power  series 
approach.  The  behavior  at  the  regular  singular  point  is  properly  accounted  for  by  using  the 
Frobenius  series  method  (Ince,  1956;  Bender,  1978)  on  the  fourth-order  equations  derived  from 
the  coupled  radial  spinor  equations.  This  is,  in  principle,  an  exact  answer  but  leaves  us  to  solve 
the  problem  of  actually  summing  the  series.  Power  series  generally  converge  slowly  and  they  are 
often  used  only  to  obtained  local  behavior,  as  in  our  previous  work.  However,  the  power  series 
solutions  can  easily  be  summed  at  large  distances  from  the  origin  if  arbitrary  precision  arithmetic 
is  used  throughout  the  calculation.  In  the  computer  algebra  programs,  Maple  and  Mathematica 
arbitrary  precision  arithmetic  is  implemented  using  the  GNU  Multiple  Precision  (GMP)  Library 
(GNU,  2010).  The  calculations  in  this  work  were  done  using  Maple  14.  Using  arbitrary 
precision,  all  of  the  terms  in  the  required  sums,  both  very  big  and  very  small,  are  accurately 
included  in  the  calculations  resulting  in  very  accurate  solutions  for  both  the  guide  wave  functions 
and  eigenvalues  that  can  be  confidently  used  in  further  detailed  investigations  of  atom  guide 
behavior. 


2.  Development  of  Radial  Equations 


The  equations  for  the  radial  wave  functions  of  the  quadrupole  guide  are  derived  by  making  use 
of  the  angular  symmetry  of  the  quadrupole  magnetic  field.  This  symmetry  is  expressed  as  the 
conservation  of  alignment,  Az  =  Lz  -Sz  (Hinds,  2001;  Lesanovsky,  2004;  Bergeman,  1989; 

Golding,  2009).  Our  approach  is  to  find  eigenfunctions  that  are  common  to  both  the 
Hamiltonian,  H  and  the  alignment  A  .  This  approach  effectively  separates  the  angular  and  radial 

dependence  of  the  problem,  leaving  only  a  coupled  system  of  radial  equations  to  solve.  These 
equations  are  solved  exactly  by  series  techniques  and  the  resulting  series  are  summed  directly 
using  arbitrary  precision  arithmetic. 

2,1  Hamiltonian  for  a  Magnetic  Atom 

The  quantum  description  of  an  atom  moving  in  a  spatially  dependent  magnetic  field  must  include 
both  internal  and  external  degrees  of  freedom.  This  means  that  the  momentum,  position,  and  spin 
degrees  of  freedom  must  all  be  treated  as  quantum  operators.  Using  m  for  the  mass  of  the  atom 

and  M  for  the  magnetic  moment,  the  Hamiltonian  is  written  as  the  sum  of  the  kinetic  energy, 

2 

A—  and  the  interaction  energy  of  the  magnetic  moment  and  the  field  B  .  The  result  is 
2m 

H  =  ^ — M-B(r).  (1) 

2m 

The  magnetic  field  5(F)  is  independent  of  the  z-direction  and  is  taken  to  be  an  ideal  quadrupole 
field  (equation  2)  extending  to  infinity  in  the  x-y  plane  and  uniform  along  z.  An  additional 
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uniform  bias  field  B0  is  added  in  the  z-direction  to  help  control  possible  spin-dependent  losses 
that  may  occur  at  the  zero  field  point  at  the  center  of  the  guide.  The  spatial  dependence  of  the 
ideal  quadrupole  field  is  given  by 


b  P)  =  Bi  {~x* + yy) + B0z  ■ 


(2) 


A  cross-sectional  plot  of  this  field  is  shown  in  figure  1 .  The  guiding  center  of  this  field  is  at  the 
origin  where  the  transverse  field  is  zero.  The  quantity  B\  is  the  magnitude  of  the  transverse  field 
gradient.  It  is  taken  as  greater  than  or  equal  to  zero  in  this  work,  although  changing  the  sign  of  B\ 
simply  changes  the  quadrupole  field  configuration  from  one  with  the  field  point  inward  along  the 
x-axis  to  one  with  the  field  pointing  inward  along  the  y- axis. 
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Figure  1 .  Cross-sectional  field  plot  of  the  general  quadrupole  guiding  field 
used  in  this  work.  Magnetic  fields  come  inwards  along  x  and  go 
outwards  along  y,  resulting  in  a  quadrupole  null  at  the  center  that 
can  be  used  for  atom  guiding. 

The  gradient  of  the  transverse  field  is  what  provides  the  trapping  force  that  keeps  the  atom 
confined  in  the  transverse  direction.  The  transverse  field  is  zero  at  the  center  of  the  guide  and  the 
only  nonzero  field  component  at  the  very  center  is  the  longitudinal  bias  field  Bo.  The  magnitude 
of  the  transverse  field  is  |f?|  =  Bl  ^ x 2  +y2  =  Btp ,  where  p  is  the  radial  coordinate  in  the 

cylindrical  system.  Even  though  the  field  is  varying  in  direction  as  one  goes  round  the  guide,  the 
contours  of  constant  field  magnitude  are  simply  circles. 
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The  potential  is  independent  of  the  position  along  the  guide.  Therefore,  pz  is  a  constant  of  the 
motion  and  is  ignored  here.  In  future  problems,  when  variations  in  the  potential  along  the  z- 
direction  are  considered,  the  variations  in  the  longitudinal  momentum  will  be  included. 

Variations  in  the  transverse  potential  should  show  up  as  variations  in  the  guide  propagation 
constant  and  this  will  affect  the  accuracy  of  sensitive  measurements  based  on  interferometry. 

In  the  model  of  the  atom  described  by  equation  1,  the  magnetic  moment  of  the  atom  is  just  the 
magnetic  moment  of  a  single  outer  electron.  The  atomic  model  used  here  can  be  thought  of  as  an 
alkali  atom  with  a  spin  zero  nucleus  and  a  total  mass  m.  The  magnetic  moment  is  M  =  yS  , 

where  y  is  the  gyromagnetic  ratio  of  the  level  considered  and  S  is  the  spin  angular  momentum  of 
the  atom.  Since  a  spin  !4  system  is  being  considered,  the  spin  angular  momentum  is  proportional 

fi 

to  the  Pauli  matrices,  S  =  —  <x. 

2 

The  Schrodinger  equation  for  the  guided  atom  eigenstate  yY ;  is 

^  £  (F)  =  (F)  (3) 


Equation  3  is  made  dimensionless  by  choosing  a  length  scale  such  that  x  —>  %x ,  where  the  new  x 

ft2  -  trk2 

is  dimensionless.  Then - V2xP  _  (r )  — > - V2lP  „  (r )  and  the  recoil  energy  is  defined  as 

trk2  1 

Er  = - ,  where  k  =  —  is  the  wavenumber  of  the  optical  field  used  to  cool  the  atom.  This  is  a 

2  M  % 

rather  arbitrary  choice  at  this  stage,  but  the  recoil  energy  is  a  convenient  scale  for  atoms  that 
have  been  laser  cooled.  After  dividing  through  by  the  recoil  energy,  pulling  out  the  relevant 


factors,  and  defining  the  dimensionless  transverse  field  parameter  bl=y 

tiB  E 

field  parameter  b0  =  y  — - ,  and  the  scaled  energy  a  =  —  . 

2Ed  Ed 


TiBx% 

~2e7 


the  longitudinal 


2.2  Eigenstates  Common  to  H  and  Az 

The  common  eigenstates  of  El  and  A_  are  found  using  the  eigenstates  of  Lz,  Sz  in  the  following 
way 


H  |  ju,e)  =  a  |  p,  a) 
Kz\p,a)  =  /u\/u,a). 


(4) 


When  the  longitudinal  bias  field,  b0  is  zero  the  Hamiltonian  is  invariant  under  various  symmetry 
operations.  This  symmetry  is  the  cause  of  a  degeneracy  of  the  eigenstates  at  b0  =  0 .  The  two 
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degenerate  states  for  each  value  of  energy  have  opposite  alignment.  The  operation  Pvox  takes  an 
eigenstate  |  ju,s)  to  its  degenerate  partner  I  -p,s)  , 

Pyax\/u,s)  =  \-/u,£),  (5) 


where  P  takes  y  to  -y  and  a  x  exchanges  the  spinor  components  when  working  in  the  position 
representation.  In  this  representation,  the  spinor  eigenfunction  is  given  by  'P  (r  )  =  (r  |  p,  s) . 
The  component  form  is 


-  JW 

(R.Ap)^2) 

y-AP’Y) 

[Rjpy"’12  J 

(6) 


and  the  degenerate  partner  is  found  as  follows 


PyPx'V  A7)  =  Py(7xeiM’ 


r+Ap) 

R-»e{p)e 


Jvll  \ 


i(pt  2 


=  e'w 


R  ApY 

R+Ap)e 


'(p!  2 


-i(p!  2 


(7) 


Once  the  eigenstate  is  found  for  positive  alignment,  the  degenerate  partner  of  negative  alignment 
is  easily  constructed.  For  this  reason,  positive  alignment  is  assumed  throughout  the  calculations 
when  b0  is  zero  and  degeneracy  is  not  an  issue. 


When  the  assumed  form  (equation  6)  is  used  in  the  Schrodinger  equation  (equation  3)  the 
angular  dependence  drops  out  leaving  only  the  following  radial  equations  for  the  specific 
alignment  // , 


,  1  (ju  +  H  2\ 

L+R+  =  d\R+  +-8pR+  ^  R+  +(e  +  b0)R+  =  bxPR, 

P  P 

LR  =82R  +-doR  -  R  +(e-b)R  =blpR+. 

P  P 


(8) 


These  two  equations  define  the  differential  operators  L+  and  L  that  will  be  used  in  section  3.  In 

the  case  of  a  spin  half  particle  with  integral  orbital  angular  momentum,  the  pair  must  be  solved 

1  3  5 

for  the  allowed  half-integral  values  of  the  alignment  //  =  +  —  ,  +  —  ,  +—•••,  although  the  solutions 
for  negative  alignment  are  simply  derived  using  symmetry  when  the  bias  field,  b0  is  zero. 
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3.  Solution  by  the  Frobenius  Technique 


The  series  solution  of  differential  equations  near  singular  points  requires  the  use  of  the  Frobenius 
method  (Ince,  1956;  Bender,  1978).  This  is  essentially  a  generalized  power  series  approach  that 
allows  an  extra  degree  of  freedom  for  proper  handling  of  the  system  behavior  near  regular 
singular  points.  Even  though  the  differential  equation  contains  singularities,  the  Frobenius 
method  produces  some  solutions  that  are  well  behaved  at  the  singular  points  and  some  solutions 
that  have  logarithmic  or  other  types  of  singular  behavior  that  do  not  usually  satisfy  the  physical 
requirements  of  a  desired  solution. 

The  differential  operator  form  of  the  coupled  radial  equations  is  useful  for  understanding  the 
series  solution  technique  that  follows.  The  operators  Z+/_  defined  in  equation  8  are  explicitly 

written  as 


P  P 


L-  =d2p+—dp-—  P  +(s-i,). 
P  P~ 


Using  operator  notation  the  coupled  system  is  just  written  as  the  pair  of  differential  equations 


L+R+=bi  PR-> 
L_R_  =blpR+. 


(10) 


3,1  Uncouple  to  Produce  Fourth  Order  Differential  Equations 

In  order  to  apply  the  Frobenius  method  to  this  system,  the  equations  are  uncoupled  producing 
two  fourth-order  equations  that  are  central  to  this  technique.  By  solving  for  both  R+  and  R 

directly  on  the  right-hand  side  in  the  above  system  (equation  10)  and  substituting  back  into  the 
other  equation,  we  find  the  pair  of  fourth-order  differential  equations 
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These  equations  are  uncoupled  in  the  sense  that  only  one  unknown  function  is  involved  in  each. 
Both  equations  are  fourth  order,  homogeneous  linear  differential  equations  and  therefore  four 
independent  solutions  are  required  for  the  general  solution  of  each.  Both  equations  have  a  regular 
singularity  at  the  origin  and  an  irregular  singularity  at  infinity.  Equations  of  this  type  are  called 
Hamburger  equations  (Ince,  1956). 

The  general  solution  of  each  equation  in  equation  1 1  is  actually  only  the  general  solution  for  one 
component  of  the  two-component  spinors  defined  in  equation  7.  The  other  component  is  found 
by  substitution  back  into  the  original  system  (equation  10).  Thus  there  are  two  distinct  ways  of 
finding  the  general  solution  of  the  system  (equation  10).  One  could  either  solve  the  first  equation 
in  equation  1 1  for  R  and  use  the  second  equation  of  equation  10  to  calculate  the  corresponding 

R  components  or  use  the  reverse  process  to  get  R  from  R  .  Because  the  Frobenius  method 

always  produces  a  log-free  solution  for  the  largest  indicial  index  found,  these  two  reversed 
solution  processes  detennine  the  two  distinct  types  of  physically  acceptable  modes  in  this 
system. 

3.2  Frobenius  Series  Approach 

The  Frobenius  series  technique  is  applied  by  assuming  a  modified  power  series  solution  of  the 
form 


Rpl=HCPnP 


n+(Jx 


(12) 


where  the  coefficients  cpn  and  the  index  cr1  are  unknown  except  that  cp0  =  1 .  In  both  the 
subscript  p\  and  in  the  coefficient  cp  ,  the  p  refers  to  the  fact  that  this  solution  process  starts 
with  the  spin-up  (or  plus)  radial  component.  The  number  1  in  the  subscripts  labels  these  solutions 
as  modes  of  type  1 .  The  reverse  solution  process  will  produce  modes  of  type  2  and  the  letter  m 
will  indicate  a  spin-down  or  (minus)  component.  Setting  R+  =  R  in  the  lower  equation  of 

equation  1 1  results  in 


f  1  A 

L  - L, 


b\P 


J 


YjCPn 


n= 0 


f  1  ^ 

L~L+ 

\  P  J 


Rpl  =blpR„  i 

oo 

p"'1  =a22>.p 


n+<r i  +1 


n= 0 


Y{cpnL_-L+-b{cpn_  1 
p 


n= 0  V 


pn+a'  =  0 


(13) 


The  result  of  letting  the  operator  L  —  L  act  on  p"+a'  is 

P 
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f  1  1 

L_-L+ 

pn+a'  = 

f  Px  P3(/I,<71) 

3  '  5 

V  P  ) 

[P  P  P  J 

p 


and  the  functions  P  are  defined  here 


(14) 


P,=z2-b02 

P3(n,<j)  =  2 eg2  +(4sn-2s -2b0)cr  +  2sn2  +(-2 b0  -2 s)n  +  s/ 2-2sp 1  +2 b0ju  +  b0 
P5(n,cr)  =  n4  +(4cr-6)n3  +(l7  /  2-lSa  —  2p2  +6cr2)n2  +  (15). 

(4a3-18o-2  +  (l7-4/r)o-  +  3/2  +  6//2  +  6//)«  + 

cr4  -6cr3  +  (-2 /T +17/2) cr2+ (3 /2  +  6//2 +  6//) cr-9//  +  //4-19/2//2-  — 


Inserting  equation  14  into  equation  13  produces 


00 

Tj(cPnPy+a'-1  +cp„P 3  {n,ax)pn^  +cp„P5  (n,at  )pn^5 -b{cpn_xpn+a' )  =  0  .  (16) 


By  adjusting  the  indices  in  the  terms  of  equation  16,  an  equation  in  which  each  term  of  the  sum 
contains  only  one  power  of  p  can  be  obtained.  In  this  form,  the  coefficients  cpn  are  defined  to  be 

zero  for  negative  n  and  it  is  easy  to  pick  out  both  the  indicial  equation  and  the  recurrence 
relations  for  the  coefficients, 

00 

YXCPnR 5  (">  ai  )  +  CPn-A  ("  "  2,  Cl  )  +  CPn_4Px  -bl2Cpn_6)p',+a'-5  =  0 

n= 0 

cpn  =  0  for  n  <  0  .  (17) 

cp0  =  l 

This  equation  is  solved  by  setting  the  factor  multiplying  p"+CT'~ 5  to  zero  for  each  n  .  This 
produces  the  following  system  of  equations  that  defines  both  the  indicial  equations  and  the 
recurrence  relations  for  the  coefficients, 

CP,A  (n,(Jx)  +  cpn_2P3  (n- 2, cr, )  +  cp^4Px  - bxcpn_6  =0,  n  =  0 . . . oo  .  (18) 


The  first  several  of  these  equations  have  been  written  out  explicitly  here, 
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n  =  0,  cp0P5  (0,cr,)  =  0, 
n  =  1,  cp,P5(l,cr,)  =  0, 

«  =  2,  cp2P5  (2,cTj)  +  c/>0P3  (0,  CTj )  =  0, 

« =  3,  c^3P5(3,cr1)  +  cp1JP3(l,cr1)  =  0, 

«  =  4,  cp4P5  (4,cTj)  +  c/>2i^  (2,cr1)  +  c/?0/J  =  0,  (19) 

n  =  5,  cp5P5  (5,  cr, )  +  cp3P3  (3,  cr, )  +  cp,F;  =  0, 
n  =  6,  cp(P5  (6, cr, )  +  cp4P2  (4,  <r, )  +  cp2Px  - b2cp0  =  0, 
n  =  7,  cp7P5  (7,cr,)  +  cp5P3  (5, cr, )  +  cp3/)  - b{cpx  =  0, 
h  =  8,  cp8Ps  (8,  crx )  +  cp6P2  (6,  ox )  +  cp4/)  -  bxcp2  =  0,  •  •  • 

Since  cp0  =  1 ,  the  n  —  0  equation  determines  the  allowed  values  of  cr,  in  terms  of  the  alignment. 
This  is  the  indicial  equation  and  is  given  byi^Cfcr,)  =  0  or 

cr,4  -6cr,3  +(-2//2  +17 / 2)cr12  +(3/2  +  6//2  +6//)cr1  -9ju  +  juA  -19/2//2  — -  =  0  (20) 
v  '  ’  16 

By  inspecting  equation  19,  it  can  be  determined  that  all  of  the  odd  n  coefficients  cpn  are  zero 
since  P5  (fcr,  I/O  and  all  odd  coefficients  are  proportional  to  cpx . 

The  roots  of  the  indicial  equation  in  equation  20  determine  the  form  of  the  solutions  for  R  x 
generated  by  the  Frobenius  technique.  The  four  roots  are 

a,  =  [p  +  l/2,-p  +  7/2,p  +  5/2,-p-l/2\ .  (21) 

The  allowed  values  for  //  are  positive  half-integers  so  all  of  the  values  of  cr,  are  integers.  In  the 

Frobenius  method,  repeated  roots  and  roots  separated  by  integers  generally  introduce  logarithmic 
singularities  into  the  solutions  and  negative  values  of  cr,  produce  solutions  with  poles  at  the 

origin.  All  of  the  singular  solutions  must  be  excluded  in  the  physical  problem.  In  the  Frobenius 
method,  the  solution  corresponding  to  the  largest  indicial  root  is  always  free  of  logarithmic 
terms.  For  //  >  1  / 2 ,  the  largest  allowed  value  of  cr,  is  always  ju  +  5/2  . 

This  calculation  is  essentially  repeated  for  the  type  2  modes  by  solving  the  first  of  equation  in 
equation  1 1  using 

oo 

*=*..=1“./’“'  <22> 

n= 0 

to  develop  the  slightly  different  indicial  equation  for  cr2 , 

a24 -6a,3  +  (-2p2  +17 /2)ct22  +(3/2  +  6p2 -6p)a2+9p  +  p4 -l9/2p2 -  —  =  0.  (23) 
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This  indicial  equation  can  be  obtained  directly  from  equation  20  by  simply  changing  the  sign  of 
//  and  the  solutions  are  easily  obtained  in  the  same  way  from  equation  2 1 , 

a2  =[-//  + 5 12, -1/2  +  p,p  +  7/2, 1/2-//].  (24) 

However,  the  largest  index  for  positive  //  is  now  //  +  7  /  2  and  since  the  largest  index  is  always 
log-free,  the  solution  (equation  22)  using  this  value  of  cr2  is  the  second  physically  acceptable 

solution  of  the  system  and  is  referred  to  as  mode  2. 

It  is  straightforward  to  see  that  for  //  >  7  /  2  there  are  at  least  two  solutions  to  equation  1 1  that  are 
unacceptable  because  of  the  pa  singularity  at  the  origin.  For  smaller  values  of  // ,  it  is  not  as 
obvious  because  there  are  some  cases  in  which  three  values  of  a  are  unique  integers.  In  these 
cases,  the  three  solutions  need  to  be  checked  more  carefully  to  be  certain  they  are  not  all  log- 
free.  We  have  checked  this  using  the  computer  algebra  program,  Maple,  and  it  turns  out  that  two, 
and  only  two,  log-free  modes  can  be  found  for  any  value  of  //  in  this  system. 


Once  the  allowed  indices  have  been  detennined,  the  coefficients  of  the  series  expansions 
equations  12  and  22  must  be  found  to  complete  the  solution.  The  cpn  are  detennined  directly 

from  equation  1 8  as 


cpn 


-cp«-2p3  p-2,fy)- cp^Px  +  b{cpn_6 
p5  l) 


n  >  2  and  even 


(25) 


and  a  similar  expression  is  obtained  by  the  same  technique  for  the  cmn .  The  recurrence  relations 
defining  the  two  physically  acceptable  solutions,  detennined  by  crl  =  p  +  5  /  2  for  the  cpn  and 
a2=  p  +  7  /  2  for  the  cmn  are 


_  bxcPn-6  +  (V  ~ g2 ) cp„_4  ~2 n(e (Ip  +  n)-b0 ) cpn_2 
n  (n  +  2) (2  p  +  n  +  3) (2  p  +  n  - 1) 

b\Cmn_6  +  (b02  - s1 ) cmn  4  -  2 (/?  + 1  +  2 //) (sn  +  £  +  b{) ) cmn_2 

CJfl  = - - - - - 

n  (4  +  (2  p  +  n  +  3)  (2  p  +  n  + 1) 

cp0  =  1,  cm0  =  1, 

cpn  =  0,  cmn  =  0  for  all  even  n  and  for  all  n  <  0 


(26) 


Using  these  recurrence  relations,  the  following  power  series  solutions  for  the  two  types  of 
physical  modes  can  be  used  directly  for  calculating  eigenvalues,  eigenfunctions,  and  the  quantum 
behavior  of  the  quadrupole  waveguide.  The  mode  1  solutions  for  both  spinor  components  are 

oo 

=  PrtV2'Zcp„p" , 

(27) 

u-1/2  co  I  p  I  U  t  jo  v  ' 

K\  =  —  Z  cPnPn  (4  F  +  2  H  n  +  6  +  5  n  +  n2 )  ■ +  — — —  p/,+3/2  Z  CP„P" 

n= 0  n= 0 
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and  the  equivalent  solutions  for  the  mode  2  spinor  components  are 


Rmi=  PM+V2YjCmnpn , 


RP2  = 


fl+ 1/2  co  (  c  —  h  \  °0 

^ ^  cmnp"  (8p  +  2p«  +  12  +  57«  +  »2)  +  - —  /U+5/2  cmnpn 

hi 


(28) 


■'l  h=0 


11=0 


4.  Series  Evaluations  Using  Arbitrary  Precision  Arithmetic 


The  waveguide  modes  are  completely  defined  by  the  coefficients  in  equation  26  and  the  series 
solutions  presented  in  equations  27  and  28.  In  order  to  find  the  eigenvalues,  the  series  solutions 
are  evaluated  at  some  large  value  of  p  and  £  is  adjusted  until  the  bound  component  of  the  spinor 
at  that  large  radius  is  near  zero.  The  bound  component  is  just  the  combination  R  +  Rm  (Hinds, 

2000;  Golding,  2009)  for  the  specific  mode  type,  either  1  or  2.  Asymptotic  analysis  shows  that 
the  bound  components  die  off  exponentially,  but  there  is  a  small  quickly  decaying  oscillatory 
component  that  can  have  an  effect  on  the  eigenvalues  calculated  using  this  technique  and  this 
possibility  must  be  studied  further. 

The  series  solutions  converge  but  they  cannot  be  easily  summed  using  the  standard  floating  point 
hardware  available  in  most  computers.  The  solution  adopted  here  is  to  use  an  existing 
implementation  of  arbitrary  precision  arithmetic  for  these  calculations.  The  computer  algebra 
program  Maple  makes  use  of  the  GNU  Multiple  Precision  (GMP)  Library  (GNU,  2010)  to 
perform  the  basic  arithmetic  calculations  addition,  subtraction,  division,  multiplication,  and 
exponentiation.  These  operations  are  all  that  are  needed  to  sum  the  power  series  solutions  in 
equations  27  and  28  to  perform  the  calculations  reported  in  section  4.1. 

4,1  Series  Calculations 

The  basic  difficulty  in  direct  summation  of  the  series  solutions  in  equations  27  and  28  is  not  that 
the  series  diverge  or  that  they  do  not  converge  quickly,  it  is  really  just  that  the  intermediate  terms 
in  the  sum  are  very  large  and  cannot  be  represented  accurately  in  hardware  floating  point. 
Arbitrary  precision  arithmetic  techniques  eliminate  this  constraint  at  the  expense  of  some  speed 
and  efficiency.  The  tradeoff  is  extremely  useful  when  extreme  accuracy  is  needed  and  other 
techniques  are  suspect. 

In  order  to  do  these  calculations  effectively,  the  number  of  digits  needed  to  represent  the  largest 
number  in  the  sum  to  at  least  the  precision  needed  in  the  final  result  should  be  estimated.  In 
addition,  the  total  number  of  terms  needed  for  convergence  at  the  desired  accuracy  needs  to  be 
known.  For  example,  in  the  solution  of  an  eigenvalue  problem,  the  value  of  an  eigenfunction 
might  be  required  to  20  significant  digits  at  p  =  40  .  In  order  to  guarantee  that  the  sums  will  be 

accurately  calculated,  individual  terms  of  the  series  should  be  calculated  for  a  reasonable  test 
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case  and  the  maximum  term  value  obtained.  If  the  maximum  value  is  lO300  and  the  desired 
precision  of  the  zero  is  1CT20  then  at  least  320  digits  are  needed  to  avoid  losing  precision  in  the 
final  sum  due  to  a  loss  of  precision  of  any  single  term  in  the  series.  In  practice,  some  guard  digits 
should  be  added  and  convergence  should  be  checked  by  changing  the  number  of  digits  and  the 
number  of  terms  used  in  the  calculation. 


Both  the  number  of  digits  and  the  number  of  terms  required  to  evaluate  a  particular  sum 
accurately  can  be  visualized  by  creating  a  term  plot.  This  is  just  a  plot  of  the  log  base  10  of  the 
magnitude  of  the  terms  in  the  sequence  to  be  summed.  The  term  plot  for  the  series  defining  Rp{ 

is  shown  in  figure  2.  The  summation  calculated  is  the  first  sum  in  equation  27.  The  terms  plotted 
are  just  the  individual  log10  \cpnp"\  for  n  from  zero  to  the  largest  value  needed.  The  number  of 

digits  needed  for  the  sum  is  just  the  difference  between  the  maximum  value  and  the  minimum 
value  assuming  the  plot  is  terminated  when  the  term  with  the  desired  precision  is  reached. 

Term  Plot  for  the  Series  Solution  of  Rp 


Figure  2.  Term  plot  of  the  log  of  the  magnitude  of  the  individual  terms 

used  in  a  calculation  of  the  first  series  in  equation  27.  The  target 
value  of  radius  is  40  and  the  maximum  eigenvalue  is  30.  The 
largest  term  in  the  sequence  is  roughly  lO100  and  therefore 
120  digits  needs  to  be  used  to  do  an  adequate  job  of  summing 
this  series  for  any  value  of  radius  up  to  40  for  eigenvalues  less  than  30. 

If  the  number  of  required  digits  is  underestimated,  the  calculation  will  breakdown  numerically  at 
values  of  the  radial  coordinate  smaller  than  the  desired  target  value.  The  results  of  using  different 
numbers  of  digits  in  the  summations  are  shown  in  figures  3,  4,  and  5.  In  figure  3,  the  number  of 
digits  used  is  15,  which  is  approximately  the  number  of  digits  used  in  floating  point  hardware 
calculations  on  many  machines.  The  calculation  breaks  down  at  around  p  =  12 .  The  range  of  the 
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calculation  can  be  extended  significantly  by  using  25  digits,  as  shown  in  figure  4.  In  figure  5, 

100  digits  were  used  in  the  calculation  and  there  are  no  obvious  problems  with  the  calculations 
over  the  whole  range  from  zero  to  40.  In  figure  6,  the  tail  of  the  100-digit  calculation  is  blown  up 
by  a  factor  of  10,000  and  a  very  clean  but  heavily  damped  oscillation  is  observed.  This  particular 
oscillation  can  be  predicted  by  asymptotic  techniques  but  it  has  been  difficult  to  verify  in 
previous  calculations  performed  at  hardware  precision.  Understanding  this  behavior  in  detail  can 
help  avoid  certain  types  of  systematic  errors  in  the  eigenvalue  calculations  that  involve 
accurately  locating  zeros  in  the  tails  of  the  eigenfunctions. 


Figure  3.  Calculation  by  series  summation  of  the  lowest  two  bound  components  using  15  digits  of 

precision.  This  kind  of  result  is  representative  of  hardware  floating  point.  Beyond  a  certain  point 
the  calculation  breaks  down.  As  the  number  of  digits  used  in  the  calculation  is  increased,  the 
target  radius  can  be  significantly  increased.  In  this  way,  hardware  floating  point  limits  the 
ability  to  obtain  accurate  eigenvalues. 
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Figure  4.  The  same  calculation  as  in  figure  3.  The  only  difference  is  that  25  digits  of  precision  are  used. 

One  can  see  that  the  calculation  can  proceed  to  much  larger  target  values  as  the  number  of  digits 
is  increased. 


Figure  5.  The  same  calculation  as  displayed  in  figures  3  and  4,  except  that  100  digits  were  used  in  this 
calculation.  The  effective  target  radius  can  be  pushed  out  to  40  with  very  clean  results.  This 
technique  lets  us  calculate  significantly  more  accurate  eigenvalues  as  well  as  perform  checks 
on  asymptotic  series  expansions  for  future  work. 
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Figure  6.  This  is  a  blown  up  version  of  the  tail  of  the  100-digit  calculation  displayed  in  figure  5.  The 
plot  is  from  radius  of  0  to  40  and  the  scale  is  expanded  by  a  factor  of  10,000.  Very  clean 
damped  oscillations  are  observed.  This  is  an  example  of  the  power  of  the  basic  arbitrary 
precision  techniques.  In  the  same  region  where  we  saw  complete  breakdown  of  the 
calculations  at  15  and  25  digits  in  figures  3  and  4,  extremely  useful  information  about  the 
asymptotic  behavior  is  obtained  by  simply  summing  the  series  using  100  digits  of 
precision.  This  cannot  be  obtained  directly  using  hardware  floating  point  techniques. 

4.2  Eigenvalues:  Shooting  or  Summing? 

In  order  to  solve  eigenvalue  problems  of  the  type  studied  here,  a  variation  of  the  shooting 
method  is  used.  The  shooting  method  is  the  technique  in  which  a  differential  equation  is 
integrated  from  a  boundary  point  where  the  solution  is  known  to  a  target  point  where  the  solution 
should  attain  a  target  value.  In  practice,  the  eigenvalue  or  some  other  parameter  can  be  adjusted 
so  that  the  unknown  function  has  the  desired  behavior  at  the  target  point.  When  the  solution  is 
acceptable,  the  eigenvalue  or  desired  value  of  the  parameter  has  been  found.  One  limitation  of 
this  technique  is  that  the  differential  equation  solver  must  be  able  to  integrate  the  system 
accurately  from  the  known  point  to  the  target  point.  This  can  be  difficult  if  the  system  is  stiff  or 
unstable  or  otherwise  very  sensitive  to  the  exact  value  of  the  eigenvalue. 

To  determine  the  eigenvalues  using  the  series  approach,  a  variation  of  the  shooting  technique  is 
used.  In  this  approach,  the  solution  does  not  depend  on  numerical  solutions  of  the  differential 
equations  from  an  initial  point  to  a  target  point.  The  only  requirement  is  to  accurately  sum  the 
series  at  the  target  point  for  various  trial  eigenvalues.  The  initial  condition  is  built  in  to  the  series 
solution.  The  eigenvalue  must  be  varied  and  the  summation  repeated  until  the  target  value  for  the 
unknown  function  is  determined.  This  type  of  effective  shooting  technique  can  be  accomplished 
over  a  large  range  of  distances  in  systems  that  support  arbitrary  precision  arithmetic  and  have 
series  solutions. 
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Because  the  sum  of  the  two  components  R  +  Rm  is  an  apparent  bound  state  for  both  mode  types 

1  and  2  (Hinds,  2000;  Golding,  2009),  the  leading  behavior  of  this  bound  component  dies  off 
exponentially  with  distance  from  the  center  of  the  guide.  Using  this  fact,  the  target  behavior  for 
the  eigenvalue  solution  is  just  that  R  +  Rm  ->  0  as  p  -»  oo  .  The  eigenvalues  are  then  readily 

found  using  a  root  solving  routine  that  uses  arbitrary  precision  arithmetic  to  look  for  zeros  of  the 
sums  defined  in  equations  27  and  28  at  large  radii.  We  have  found  that  in  this  way  the  target 
point  can  be  much  farther  from  the  origin  than  it  can  using  the  similar  ODE  based  shooting 
techniques.  This  results  in  much  more  precise  calculations  of  the  eigenvalues.  This  enhanced 
resolution  is  expected  to  be  useful  as  these  functions  are  used  for  further  calculations  on  the 
detailed  properties  of  atom  guides. 

The  results  of  eigenvalue  calculations  for  //  =  1  /  2,  b0  =  0  and  bl  =  1  are  shown  in  figure  7.  The 

first  52  eigenvalues  are  shown.  The  eigenvalue  data  can  be  approximated  by  the  simple  equation 
1.893790658  +  1. 270 128257/? 0  7355 .  In  table  1,  the  first  26  values  of  each  type  of  mode  are 

presented.  Careful  inspection  shows  that  eigenvalues  form  an  interleaved  series  in  the  sense  that 
the  first  eigenvalue  is  a  type  1  mode  and  second  is  a  type  2  mode;  this  pattern  holds  up  as  far  as 
we  have  checked.  In  figure  7,  mode  type  1  points  are  red  and  mode  type  2  points  are  blue. 
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Figure  7.  Plot  of  the  first  52  eigenvalues.  The  red  points  are  mode  type  1  eigenvalues  and  the  blue  points  are 
mode  type  2  eigenvalues;  an  approximate  fit  to  the  eigenvalues  as  a  function  of  n  is  given  in  the 
text. 
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Table  1.  First  52  energy  eigenvalues  for  our  system  with  alignment  1/2, 

zero  bias  field  and  a  transverse  gradient  ofbl=l.  The  sequence  of 
eigenvalues  alternates  between  type  1  and  type  2  eigenvalues  as  n 
is  increased.  Because  of  this  double  counting  the  index  n  in  the 
table  is  not  the  same  as  the  index  n  used  in  the  plot  of  figure  7. 


n 

Type  1  Eigenvalue 

Type  2  Eigenvalue 

i 

1.893790658 

2.770959115 

2 

3.735223647 

4.436734695 

3 

5.211244681 

5.827843267 

4 

6.504703191 

7.067475471 

5 

7.681587325 

8.20585412 

6 

8.774989674 

9.269728735 

7 

9.80443751 

10.27548299 

8 

10.78265521 

11.23407047 

9 

11.71851157 

12.15327128 

10 

12.61849734 

13.03886265 

11 

13.48754237 

13.89528345 

12 

14.32950179 

14.72603767 

13 

15.14746211 

15.53395267 

14 

15.9439429 

16.32135181 

15 

16.72103469 

17.09017376 

16 

17.48049623 

17.84205754 

17 

18.22382493 

18.57840462 

18 

18.952309 

19.30042528 

19 

19.66706684 

20.00917379 

20 

20.36907731 

20.70557569 

21 

21.05920333 

21.39044912 

22 

21.73821055 

22.06452178 

23 

22.40678229 

22.72844459 

24 

23.06553163 

23.38280277 

25 

23.7150113 

24.02812491 

26 

24.35572179 

24.66489048 

4.3  Eigenfunctions 

Once  the  eigenvalues  are  determined,  the  corresponding  eigenfunctions  are  found  by  a 
straightforward  summation  of  the  series  solutions  equations  27  and  28  at  any  desired  value  of  p  . 

The  eigenfunctions  determined  in  this  way  display  detailed  behavior  in  regions  beyond  those 
accessible  using  the  hardware  precision  ODE  shooting  techniques.  This  type  of  behavior  is 
shown  in  figures  8  and  9.  It  is  clear  that  the  bound  component  does  not  simply  decay  away 
exponentially  but  has  a  small  and  quickly  decaying  oscillatory  component  as  well.  This  behavior 
is  explained  by  a  detailed  asymptotic  solution  of  the  problem.  The  ability  to  compare  the  direct 
calculations  of  the  series  technique  with  the  asymptotic  solutions  over  wide  ranges  will  be  very 
helpful  in  checking  various  aspects  of  the  asymptotic  behavior  of  these  solutions. 
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Figure  8.  The  first  10  type  1  mode  radial  eigenfunctions.  The  black  curve  with  no  zero  crossings  is  the  ground  state 
wave  function.  The  first  excited  state  is  the  blue  curve. 


Figure  9.  The  first  10  type  2  mode  radial  eigenfunctions.  The  black  curve  is  the  lowest  energy  type  2  mode  but  is 
actually  the  first  excited  state  of  the  system. 
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5.  Conclusions 


The  use  of  series  solution  techniques  along  with  arbitrary  precision  arithmetic  is  a  very  powerful 
combination  for  investigating  the  solutions  of  differential  equations  in  regions  where  ordinary 
hardware  precision  numerical  integration  techniques  break  down.  The  use  of  arbitrary  precision 
in  differential  equation  solvers  is  implemented  in  Maple  and  these  can  be  made  to  produce 
accurate  solutions  at  large  radii.  However,  the  solution  of  a  differential  equation  by  the  series 
technique  simply  requires  the  appropriate  power  series  formulation  of  the  solution  and  a  software 
system  that  implements  an  arbitrary  precision  library,  such  as  GMP  (GNU,  2010).  By  using  a 
judicious  choice  of  the  number  of  digits  and  terms  needed  for  a  particular  set  of  calculations,  the 
solutions  can  be  calculated  very  quickly  and  with  excellent  resolution.  We  have  found  this 
technique  to  be  extremely  useful  for  calculating  radial  wavefunctions  and  eigenfunctions  for  the 
atomic  waveguide  problem  as  it  produces  accurate  solutions  in  situations  where  spinor  equations 
can  be  uncoupled  to  produce  higher  order  differential  equations  that  are  difficult  to  solve 
analytically. 
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